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Pore  network  simulations  are  performed  to  study  water  transport  in  gas  diffusion  layers  (GDLs )  of  polymer 
electrolyte  membrane  fuel  cells  (PEMFCs).  The  transport  and  equilibrium  properties  are  shown  to  be 
scale  dependent  in  a  thin  system  like  a  GDL.  A  distinguishing  feature  of  such  a  thin  system  is  the  lack  of 
length  scale  separation  between  the  system  size  and  the  size  of  the  representative  elementary  volume 
(REV)  over  which  are  supposed  to  be  defined  the  macroscopic  properties  within  the  framework  of  the 
continuum  approach  to  porous  media.  Owing  to  the  lack  of  length  scale  separation,  two-phase  flow 
traditional  continuum  models  are  expected  to  offer  poor  predictions  of  water  distribution  in  a  GDL.  This 
is  illustrated  through  comparisons  with  results  from  the  pore  network  model.  The  influence  of  inlet 
boundary  conditions  on  invasion  patterns  is  studied  and  shown  to  affect  greatly  the  saturation  profiles. 
The  effects  of  GDL  differential  compression  and  partial  coverage  of  outlet  surface  are  also  investigated. 

©  2009  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Fuel  cells  are  attracting  increasing  attention  as  an  important 
energy  converter  and  many  groups  work  throughout  the  world 
in  order  to  improve  their  performances,  efficiency,  reliability, 
manufacturability  and  cost-effectiveness.  In  particular,  polymer 
electrolyte  membrane  (PEM)  fuel  cells  are  under  intense  studies 
and  development.  One  well  identified  limitation  of  the  PEMFC  per¬ 
formances  at  high  current  densities  originates  in  the  blockage  of 
the  pores  of  the  GDL  (Gas  diffusion  layer)  by  liquid  water.  This 
blockage  constraints  the  reactant  transport  to  the  catalyst  layers. 
This  problem  has  motivated  many  studies  on  the  modelling  of 
two-phase  transport  in  the  GDL.  One  can  distinguish  three  main 
approaches:  continuum  models,  pore  network  models,  direct  sim¬ 
ulations.  Continuum  models  refer  to  the  traditional  models  widely 
used  in  many  applications  involving  porous  media;  these  notably 
include  petroleum  engineering  and  subsurface  hydrology.  The 
porous  medium  is  treated  as  a  hypothetical  effective  continuum  and 
the  models  involve  volume-averaged  quantities  such  as  saturation 
and  rely  on  phenomenological  relationships  such  as  the  generalized 
Darcy’s  law.  Key  parameters  in  this  approach  are  the  permeability, 
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relative  permeabilities  and  the  capillary  pressure-saturation  curve. 
Despite  recurrent  questions  about  the  relevant  functional  forms  to 
be  used  for  these  parameters  for  the  fibrous  media  forming  the 
GDLs,  this  approach  has  been  used  extensively  for  modelling  water- 
gas  flows  in  GDLs,  e.g.  [1-3]  and  references  therein.  As  discussed 
for  example  in  [4],  there  are  however  serious  concerns  about  using 
this  type  of  model  for  a  thin  system  like  a  GDL,  especially  due  to 
the  lack  of  a  significant  length  scale  separation  between  the  pore 
scale  and  the  GDL  thickness.  This  has  motivated  in  part  the  two  main 
alternative  approaches  listed  above.  Direct  simulations  aim  at  com¬ 
puting  the  two-phase  flow  by  solving  directly  the  problem  at  the 
pore  scale.  This  can  be  performed  using  for  example  a  Lattice  Boltz¬ 
mann  model,  [5].  Although  interesting  for  a  better  understanding 
of  two-phase  flow  in  a  complex  microstructure  such  as  the  one  of  a 
GDL,  direct  simulations  are  of  limited  practical  use  because  of  the 
large  computing  resources  needed.  In  this  respect,  pore  network 
models  are  much  more  efficient  and  appear  as  a  good  option  not 
only  for  determining  the  parameters  needed  in  the  traditional  con¬ 
tinuum  approach,  e.g.  [6],  but  also  as  a  valuable  tool  for  gaining  a 
better  understanding  of  the  phenomena  occurring  at  the  pore  net¬ 
work  scale  as  well  as  for  possibly  improving  the  design  of  GDLs,  e.g. 
[  7-9  ]  and  references  therein.  In  this  paper,  we  also  use  pore  network 
models.  The  main  objective  is  to  evaluate  what  can  be  expected 
from  continuum  models  and  also  to  study  several  effects  on  water 
transport  in  GDL.  These  include  the  scale  effect  associated  with  the 
thin  nature  of  GDL,  the  injection  condition  at  the  GDL  inlet  and  the 


0378-7753 /$  -  see  front  matter  ©  2009  Elsevier  B.V.  All  rights  reserved. 
doi:10.1016/j.jpowsour.2009.02.090 


M.  Rebai,  M.  Prat  /  Journal  of  Power  Sources  192  (2009)  534-543 


535 


Fig.  1.  (a)  Sketch  of  GDL  unit  cell,  (b)  The  two  injection  conditions  considered  in  the  paper:  GDL  in  contact  with  a  reservoir  at  uniform  pressure,  injection  from  one  single 
pore  located  in  the  centre  of  GDL  inlet  surface.  These  boundary  conditions  are  referred  to  as  the  surface  injection  b.c.  and  one  point  injection  b.c.,  respectively. 


screen  effect  due  to  the  fact  that  only  a  fraction  of  the  GDL  surface 
adjacent  to  the  bi-polar  plate  is  in  contact  with  the  gas  channels. 
Impact  of  possible  effects  of  GDL  differential  compression  is  also 
investigated. 

The  paper  is  organized  as  follows:  first  we  recall  in  Section  2 
some  basic  results  on  drainage  in  porous  media  so  as  to  determine 
the  most  likely  two-phase  flow  regime  in  a  GDL.  A  brief  discus¬ 
sion  about  water  injection  in  the  GDL  is  presented  in  Section  3.  In 


Fig.  2.  Schematic  of  cubic  pore  network.  The  size  of  the  network  shown  here  is 
6x6x2  whereas  a  40  x  40  x  N  with  N  e  [4,40]  is  used  for  the  simulations  presented 
in  the  present  study.  For  simplicity,  the  pore  size  and  the  throat  size  shown  here  are 
uniform.  Actually  the  throat  and  pore  sizes  are  randomly  distributed  as  explained  in 
the  text. 


Section  4,  we  concentrate  on  finite  size  effects  in  a  system  of  large 
aspect  ratio  presumably  representative  of  a  GDL  unit  cell.  In  Section 
5,  we  explore  the  impact  of  screen  and  compression  effects.  Section 
6  is  devoted  to  a  brief  investigation  of  a  different  scenario  in  which 
water  enters  the  GDL  from  an  isolated  source  rather  than  from  all 
the  pores  in  contact  with  the  MPL  (Microporous  layer,  which  is  con¬ 
sidered  as  a  layer  distinct  from  the  GDL  in  this  paper  and  not  as  a  part 
of  the  GDL)  or  the  Catalyst  layer  (CL).  Some  comparisons  between 
the  continuum  model  and  pore  network  simulations  are  presented 
in  Section  7.  Implications  of  the  results  as  well  as  possible  impacts  of 
additional  effects  (mixed  wettability,  GDL  anisotropy)  are  discussed 
in  Section  8. 

Before  going  into  the  details  of  the  various  sections,  the  configu¬ 
ration  and  the  assumptions  considered  are  presented.  Throughout 
this  paper,  we  assume  for  simplicity  that  the  porous  matrix  is 
hydrophobic  (as  in  many  previous  studies)  and  that  the  GDL 
microstructure  is  isotropic  (assumption  often  taken  with  the  con¬ 
tinuum  models).  We  are  aware  that  these  assumptions  can  be 
questioned  and  this  will  be  briefly  considered  in  Section  8. 

Also  we  do  not  consider  all  the  GDL  but  restrict  our  attention  to 
a  sub-region  of  the  GDL  referred  to  as  a  unit  cell  or  a  representative 
domain,  see  Fig.  1 .  Noting  that  the  GDL  is  in  contact  with  the  bipolar 
plate  in  which,  at  least  for  the  most  common  designs,  regularly 
spaced  channels  are  etched,  we  take  the  distance  L  between  two 
channels  as  representative  size  of  a  GDL  unit  cell  in  the  in-plane 
direction.  Hence  as  a  reasonable  unit  cell  of  the  system  GDL/bipolar 
plate,  we  consider  a  porous  domain  of  size  Lx  Lx  l  where  L  is  the 
GDL  thickness.  Representative  values  of  L  and  l  are  L  ^  2  mm  and 
i  «  [170-400]  |jim,  e.g.  [10].  The  diameter  of  fibres  forming  the  GDL 
is  typically  of  the  order  of  10  p,m  whereas  the  mean  pore  size  is  on 
the  order  of  50-60  p,m,  [6,11  ].  Hence,  measured  in  pore  size,  a  GDL 
is  typically  less  than  10  pore  sizes  thick  whereas  L  is  on  the  order 
of  40  pore  sizes.  Compared  to  most  other  porous  materials,  a  very 
particular  feature  of  GDL  is  thus  to  be  a  thin  system,  not  so  much 
in  terms  of  aspect  ratio,  i.e.  i  <  L,  but  because  there  are  only  a  few 
pores  across  the  GDL  thickness. 

As  discussed  in  [4],  the  details  on  how  water  is  produced  within 
the  catalyst  layer  and  reaches  the  GDL  are  still  somewhat  unclear. 
As  in  most  of  the  previous  works  on  the  subject,  we  assume  that 
water  enters  the  GDL  in  liquid  phase  from  the  surface  of  the  GDL 
in  contact  with  either  the  microporous  layer,  if  any,  or  the  catalyst 
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layer.  We  assume  that  the  GDL  is  dry  initially  and  we  are  mainly 
interested  in  the  water  invasion  process  in  this  initially  dry  GDL 

As  mentioned  before  the  study  is  based  on  pore  network  simu¬ 
lations.  For  simplicity,  we  consider  a  cubic  network.  As  sketched  in 
Fig.  2,  pores  of  cubic  shape  are  regularly  placed  on  a  3D  cartesian 
grid  (with  a  denoting  the  lattice  spacing).  Two  first  neighbour  pores 
are  linked  by  a  channel  of  square  cross-section.  Such  a  channel  is 
referred  to  as  a  bond  or  throat.  The  pore  size  dp  (respectively  the 
throat  cross  section  side  length  dt),  which  corresponds  to  the  diam¬ 
eter  of  the  largest  sphere  inscribed  within  the  pore,  is  randomly 
distributed  according  to  a  uniform  probability  law  in  the  range 
[^pmin*  dpmax]  (respectively  in  the  range  [dtmin,  dtmax]  with  the  con¬ 
straint  dtmax  <  dpmin).  Simulations  presented  in  what  follows  have 
been  performed  with  a  =  50[xm,  dpmax  =  46.5  [xm,  dpmin  =  34  |jim, 
dtmax  =  34  pan,  dtmin  =  20  pan.  These  values  are  representative  of 
pore  size  distributions  in  GDL,  e.g.  [6]. 

The  network  size  is  characterized  by  the  number  of  pores  placed 
in  each  direction.  In  accordance  with  the  representative  values 
given  above,  we  consider  40  x  40  x  N  networks  with  N  varying  in 
the  range  [4,40]  (with  a  lattice  spacing  of  50  pan,  N= 4  corresponds 
for  example  to  a  GDL  thickness  of  200  pm). 

The  water  invasion  process  considered  in  what  follows  (see  Sec¬ 
tion  2)  is  sensitive  to  the  particular  realization  of  the  network 
considered  (since  the  pore  and  throat  sizes  are  random  variables,  a 
realization  of  the  network  refers  to  one  particular  random  drawing 
of  pore  and  throat  sizes).  As  a  result,  ensemble  averaging  is  needed 
to  obtain  the  mean  behaviour  of  a  given  variable  of  interest.  Unless 
otherwise  mentioned,  100  such  realizations  of  each  considered  net¬ 
work  have  been  typically  generated  to  obtain  the  results  presented 
in  this  paper. 


2.  Theory  of  drainage  and  thin  systems 


When  water  displaces  a  gas  in  a  hydrophobic  system,  water  is 
the  non-wetting  fluid  and  the  gas  is  the  wetting  one.  In  the  porous 
media  literature,  e.g.  [12],  such  a  process,  i.e.  the  immiscible  dis¬ 
placement  of  a  wetting  fluid  by  a  non-wetting  one  is  referred  to  as 
drainage.  When  the  gravity  effects  can  be  neglected,  a  reasonable 
assumption  in  the  present  context,  three  main  asymptotic  invasion 
regimes  can  be  distinguished  depending  on  the  values  of  the  capil¬ 
lary  number  Ca  (which  is  the  ratio  between  viscous  forces  acting  at 
the  pore  scale  in  water  and  capillary  forces,  Ca  =  /xnwU/y  cos  0W, 
U  is  a  reference  velocity,  where  y  is  the  surface  tension,  /xnw  is 
the  liquid  water  dynamic  viscosity,  0W  is  the  contact  angle  taken 
in  the  wetting  fluid  (air))  and  the  ratio  of  dynamic  viscosities 
of  the  two  fluids  M,  M=/iwl/i nw~10-3  for  the  water/air  system, 
[13,14].  These  regimes  are  the  invasion  percolation  (IP)  regime 
(very  low  Ca),  the  capillary-viscous  fingering  regime  (very  high 
M  and  Ca)  and  the  compact  regime  (very  low  M  and  high  Ca).  As 
discussed  in  [4],  the  capillary  number  can  be  expressed  here  as 
Ca  =  finwIMH20/y  cosOw2Fpnw  where  F  is  the  Faraday’s  constant, 
Mh2o  is  the  water  molecular  weight,  pnw  is  the  liquid  water  den¬ 
sity  and  I  is  the  current  density.  With  /=  1 A  cm-2  and  0W  =  70°  (see 
[4]  for  a  discussion  about  the  value  of  contact  angle),  this  leads  to 
Ca  ~  10-7  to  10-8,  a  low  value  presumably  consistent  with  the  IP 
regime.  However  as  pointed  out  in  [14]  and  [15],  the  estimate  of 
the  capillary  number  is  not  sufficient  to  determine  which  regime  is 
prevailing.  As  discussed  in  [14],  the  IP  regime  is  expected  only  for 
system  of  size  lower  than  the  length  Xe  (measured  in  lattice  unit  a) 
which  is  given  by, 


c^?X1+^+v^D_1^v  ~  ft 
E 


(1) 


where  D  =  2.52  (D  is  the  mass  fractal  dimension  of  the  percolation 
cluster),  v  =  0.88,  £  =  1.12  (v  and  £  are  the  conductance  and  corre¬ 
lation  exponents  of  the  percolation  theory,  [16]),  /3  «  0.01,  E  is  the 


Fig.  3.  Evolution  of  capillary  number  Ca  marking  the  transition  between  a  pure  IP 
regime  and  the  regime  IPSG  where  viscous  effects  cannot  be  ignored  as  a  function 
of  pore-network  thickness. 


dimensionless  standard  deviation  of  the  throat  size  distribution.  For 
a  uniform  distribution,  S  =  2(dtmax  -  dtmin)/Vl2(dtmax  +  dtmin). 
In  Eq.  (1),  c  is  a  constant  which  has  be  computed  for  our  network 
and  is  found  to  be  c^40.  As  shown  in  [14]  or  [15],  Eq.  (1 )  is  obtained 
by  estimating  the  evolution  of  capillary  pressure  induced  by  vis¬ 
cous  effects  along  an  IP  cluster  and  corresponds  here  to  the  special 
case  where  the  pressure  drop  in  the  less  viscous  fluid  (the  gas)  is 
neglected.  From  Eq.  (1),  we  can  estimate  the  evolution  of  capillary 
number  above  which  the  displacement  ceases  to  be  a  pure  IP  dis¬ 
placement.  The  results  are  shown  in  Fig.  3  where  IPSG  stands  for 
invasion  percolation  in  a  stabilizing  gradient  and  corresponds  to 
the  regime  where  viscous  effects  are  non-negligible.  As  can  be  seen 
from  Fig.  3,  the  values  of  Ca  10-7  to  10-8  estimated  before  corre¬ 
spond  in  fact  to  the  values  marking  the  limit  of  the  IP  regime  for 
N<  10,  which  is  the  range  of  N  expected  for  the  GDL.  Hence  this 
confirms  that  the  IP  model  is  a  valuable  tool  to  simulate  the  dis¬ 
placement  of  gas  by  water  in  a  GDL.  However,  one  should  keep  in 
mind  that  the  system  operates  close  to  the  limit  of  validity  of  the 
IP  model.  Therefore,  viscous  effects  might  become  non-negligible 
if  higher  liquid  flow  rates  were  to  be  considered  (because  of,  for 
instance,  improvements  in  the  fuel  cell  performances). 

3.  Injection  conditions 

As  mentioned  before,  the  exact  conditions  under  which  water 
is  produced  in  the  active  layer  and  reaches  the  GDL  (after  travel¬ 
ling  across  the  microporous  layer  (MPL)  if  any)  are  still  somewhat 
unclear,  e.g.  [1,17].  As  in  most  previous  investigations,  we  assume 
here  that  water  enters  the  GDL  in  liquid  phase  and  that  condensa¬ 
tion  of  water  vapour  can  be  neglected.  In  the  traditional  application 
of  IP  model,  the  entrance  surface  of  the  porous  medium  is  supposed 
to  be  in  contact  with  a  reservoir  of  non-wetting  fluid  at  uniform 
pressure.  However,  the  situation  may  be  different  for  a  GDL  in  an 
operating  fuel  cell  owing  to  the  presence  of  the  adjacent  porous 
layer,  i.e.  the  active  layer  or  the  MPL.  As  interestingly  pointed  out  in 
[9]  and  [17],  more  realistic  conditions  could  be  multiple  injections 
through  each  pore  of  the  GDL  entrance  surface  or  only  through  some 
of  them.  Since  the  traditional  condition  is  widely  used  and  is  pre¬ 
sumably  representative  of,  at  least,  ex-situ  experiments,  e.g.  [18], 
most  of  our  simulations  will  rely  on  this  condition.  However,  we 
will  consider  in  Section  6  also  the  extreme  case  of  the  injection 
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through  one  single  injection  pore  so  as  to  have  some  indication  of 
the  possible  impact  of  boundary  condition  on  invasion  pattern  and 
pore  occupancy  at  breakthrough.  These  two  cases  are  sketched  in 
Fig.  1  and  are  referred  to  as  the  “surface  injection”  condition  and 
the  “one  point  injection  “condition,  respectively. 

4.  Finite  size  effects  in  thin  system 

One  major  feature  of  one  invasion  percolation  pattern  is  to  be  a 
fractal  pattern.  As  a  result,  many  quantities  are  scale  dependent, 
especially  at  breakthrough  (BT),  i.e.  when  water  reaches  for  the 
first  time  the  GDL/bipolar  plate  interface.  This  scale  dependence 
is  illustrated  here  through  the  consideration  of  several  quantities 
of  interest,  ignoring  throughout  this  section  the  screen  effect  men¬ 
tioned  in  the  introduction,  i.e.  the  water  injected  on  one  side  of  the 
GDL  is  free  to  exit  from  any  throat  present  on  the  GDL  opposite  sur¬ 
face  (with  the  notation  used  in  Section  5,  this  corresponds  to  0  =  1, 
where  0  is  the  fraction  of  the  GDL  unit  cell  surface  on  the  bipolar 
plate  side  in  contact  with  the  gas  channel). 

4  A.  Overall  saturation  at  breakthrough 

The  overall  water  saturation  SBT  at  breakthrough  is  the  volume 
fraction  of  pore  space  in  the  40  x  40  x  N  network  which  is  occupied 
by  water  at  breakthrough  whereas  (S)BT  is  the  overall  water  sat¬ 
uration  at  breakthrough  averaged  over  many  realizations  (500  for 
N = 4, 350  for  N  =  10  and  250  for  N  >  10).  To  compute  SBT  for  a  particu¬ 
lar  realization,  we  use  the  IP  algorithm  without  trapping  (trapping 
is  negligible  in  3D  network  up  to  breakthrough),  as  described  for 
example  in  [4].  For  a  porous  domain  LxLxL  where  L  is  measured 
in  lattice  unit,  (S)Bt  scales  as,  e.g.  [16], 

(5)bT  0(1  £0.48 

Hence  the  fraction  of  pore  space  1  -  (S)Bt  left  available  for  the 
gaseous  reactant  to  reach  the  catalyst  layer  increases  with  the  size 
of  the  domain  (measured  in  lattice  spacing  unit). 

Consider  now  a  LxLxN  pore  network  of  fixed  lateral  extension 
L  (1  =  40  in  our  case)  and  study  the  influence  of  N  on  (S)Bt.  Pore 
network  simulations  based  on  the  IP  algorithm  lead  to  the  results 
shown  in  Fig.  4  (“surface  injection”  results  in  Fig.  4).  As  can  be  seen, 
(S)Bt  increases  with  N.  Interestingly,  it  can  be  seen  from  Fig.  4  that 
(S)BT  is  very  low,  close  to  4%,  for  the  thinnest  system  considered  (4 
pore  sizes  thick)  and  increases  rapidly  with  thickness  up  to  13%  for 
the  40  pore  sizes  thick  system  (for  which  (S)BT  is  given  by  Eq.  (2) 
with  L  =  40).  Hence  these  results  indicate  that  the  thinnest  the  GDL 
(expressed  in  lattice  unit),  the  smallest  the  fraction  of  blocked  pores 
at  breakthrough.  Interestingly,  the  fraction  of  pore  volume  available 
to  the  gas  can  be  considered  as  high  whatever  the  GDL  thickness  in 
the  range  of  thickness  investigated,  which  is  a  clear  indication  that 
the  IP  regime  is  well  adapted  for  draining  the  water  across  the  GDL 
with  a  high  degree  of  accessibility  to  the  catalyst  layer  for  the  gas. 
Also,  a  typical  feature  of  thin  system  is  the  significance  of  statistical 
fluctuations  as  exemplified  by  the  error  bars  in  Fig.  4. 

4.2.  Capillary  pressure  curve 

As  mentioned  in  the  introduction,  the  capillary-pressure  curve 
Pc(5)  (or  equivalently  Pc(Sg)  where  5g  =  l  -  5)  is  one  of  the  key 
parameters  entering  into  the  traditional  continuum  models.  As 
exemplified  for  example  in  [6],  pore  network  models  can  be  used 
to  determine  PC(S).  To  this  end,  we  used  the  algorithm  described 
hereafter,  which  also  relies  on  IP  concepts.  Due  to  the  random  fluc¬ 
tuations,  Pc(5)  has  to  be  determined  over  many  realizations  of  the 
network.  For  a  given  realization,  we  start  with  a  dry  network  and 
determine  the  overall  capillary  pressure  as  a  function  of  S  for  suc- 


Fig.  4.  Evolution  of  mean  overall  water  saturation  at  breakthrough  as  a  function 
of  network  thickness  for  the  two  boundary  conditions  depicted  in  Fig.  1.  The  error 
bars  represent  ±1  standard  deviation  around  the  mean  values  over  the  number  of 
realizations  considered. 

cessive  states  of  hydrostatic  equilibrium  corresponding  to  small 
increment  dPnw  in  the  non-wetting  fluid  pressure  (the  wetting  fluid 
pressure  is  maintained  constant).  The  algorithm  used  to  determine 
the  saturation  evolution  right  after  a  non-wetting  fluid  pressure 
increment  reads: 

1.  Identification  of  each  bond  that  can  be  invaded,  that  is  each  wet¬ 
ting  fluid  bond  such  that  its  capillary  pressure  threshold  pc  is 
lower  than  the  pressure  difference  Pnw  -  Pw  between  the  two 
fluids. 

2.  Identification  of  clusters  formed  by  the  bonds  identified  in  (1) 
using  the  Hoshen-Kopelman  algorithm  [19]. 

3.  Identification  of  clusters  among  the  clusters  identified  in  (2)  in 
contact  with  the  invading  non-wetting  phase. 

4.  Invasion  of  all  clusters  identified  in  (3). 

5.  Computation  of  overall  saturation. 

Note  that  trapping  is  not  taken  into  account  here  as  well  as  in 
the  rest  of  the  paper,  see  for  example  [20]  and  references  therein  for 
more  details  on  this  aspect.  The  local  capillary  pressure  threshold 
pc  of  a  throat  is  expressed  using  the  Young-Laplace  equation  as 

pc.i>^L  (3) 

at 

Using  this  algorithm  for  many  realizations  ( see  Section  4.1  for  the 
details  on  the  number  of  realizations  generated)  leads  to  the  results 
shown  in  Fig.  5.  As  can  be  seen  from  Fig.  5,  the  capillary  pressure 
curve  is  thickness  dependent  but  this  dependence  is  only  marked 
in  the  region  of  high  gas  saturation,  i.e.  near  breakthrough.  The  evo¬ 
lution  of  (Pc(Sg))  in  this  range  of  saturation  tends  to  be  steeper  as 
N  increases.  A  detailed  view  of  this  evolution  is  shown  in  Fig.  5b 
together  with  error  bars  corresponding  to  ±1  standard  deviation 
around  the  mean  values.  These  results  suggest  that  the  experimen¬ 
tal  measurement  of  capillary  pressure  curve  in  a  thin  system  should 
be  performed  with  the  thin  system  as  it  is,  i.e.  determining  for  exam¬ 
ple  Pc(Sg)  for  a  thick  system  formed  by  a  sandwich  of  several  layers 
of  the  thin  system  is  not  a  good  idea  if  an  accurate  description  of 
Pc(5g)  for  Sg  close  to  1  is  sought.  This  region  of  the  capillary  curve  is 
expected  to  be  a  very  sensitive  region  for  the  continuum  model  if 
one  wishes  to  describe  the  invasion  of  water  in  a  dry  hydrophobic 
GDL  with  this  type  of  model  (this  is  not,  however,  the  main  prob- 
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Fig.  5.  (a)  Influence  of  network  thickness  on  capillary  pressure  curve,  (b)  Detailed  view  of  capillary  pressure  curve  for  the  high  overall  gas  saturation,  the  error  bars  represent  ±1 
standard  deviation  around  the  mean  values  over  the  number  of  realizations  considered.  Pcref  is  the  capillary  pressure  threshold  given  by  Eq.  (3)  for  dt  =  dt  =  (dtmax  +  dtmin/2). 


lem  with  the  continuum  models  as  discussed  in  Section  7).  This  is 
due  to  the  fact  that  the  derivative  dPc/dS  is  one  of  the  parameters 
of  the  continuum  models.  As  can  be  seen  from  Fig.  5,  dPc/dS  tends 
to  diverge  as  Sg  ->  1.  This  issue  is  however  beyond  the  scope  of  the 
present  paper  and  will  be  addressed  in  a  forthcoming  work  (see 
however  Sections  7  and  8  below). 

4.3.  Permeability  and  relative  permeabilities 

The  intrinsic  permeability  and  the  relative  permeabilities  are 
also  important  parameters  of  continuum  models.  To  determine 
the  permeability  and  relative  permeabilities  of  the  network,  local 
hydraulic  conductances  must  be  assigned  to  each  throat  and  pore. 
As  in  [6],  the  local  conductance  of  a  throat  or  a  pore  is  expressed  as 

d4 

8~  14.03 (  (  ^ 

where  d  is  the  pore  or  throat  size  and  l  is  the  pore  or  throat  length. 
The  local  flow  rate  q  between  two  pores  is  expressed  as 

q  =  —  AP  (5) 

fi 

where  pt  is  the  dynamic  viscosity  of  the  considered  fluid  and  gT  is 
the  hydraulic  conductance  between  the  two  pores  obtained  as  the 
harmonic  average  of  the  conductances  of  the  two  half  pores  and  the 


throat  over  which  there  exists  the  pressure  difference  AP.  The  pro¬ 
cedure  is  exactly  similar  to  the  one  described  in  [6],  where  more 
details  can  be  found.  Expressing  the  mass  conservation  equation 
over  each  pore  yields  a  system  of  linear  equations  that  is  solved 
numerically.  As  boundary  conditions,  arbitrary  pressures  Pen  and 
Pex  are  prescribed  on  inlet  and  outlet  surfaces  of  network.  This 
yields  the  total  flow  Q  over  the  network.  Once  Q.  is  determined, 
the  permeability  I<  of  the  network  is  found  from  Darcy’s  law: 


K  Pen  —  Pex 
pt  i 


(6) 


In  our  case,  this  gives  (K)  ^4.3  x  10-11  m2  without  significant 
variations  with  the  size  N  of  network. 

For  determining  the  relative  permeabilities  (/<r),  we  combine  the 
procedure  described  in  Section  4.2  to  determine  the  fluids  distri¬ 
bution  within  the  network  and  the  computation  of  permeability 
described  above.  For  example  the  effective  permeability  I<  kr  of  one 
fluid  is  determined  using  the  same  procedure  as  for  the  intrinsic 
permeability  except  that  only  the  sub-network  corresponding  to 
this  fluid  is  considered.  Again  more  details  on  this  type  of  compu¬ 
tation  can  be  found  in  [6]. 

Not  surprisingly,  Fig.  6  shows  that  the  relative  permeabilities 
(averaged  over  the  same  number  of  realizations  as  for  (S)bt  and 
(Pc(Sg)))  are  network  thickness  dependent.  Again  the  change  is 
quite  significant  when  the  network  size  varies  from  4  to  10  and 


Fig.  6.  (a)  Influence  of  network  thickness  on  relative  permeability  curves,  (b)  Detailed  view  for  the  high  overall  gas  saturation,  the  error  bars  represent  ±1  standard  deviation 
around  the  mean  values  over  the  number  of  realizations  considered. 
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Fig.  7.  Definition  of  various  zones  for  the  analysis  of  screen  and  differential  com¬ 
pression  effects.  Zones  1  and  2  can  be  affected  by  mechanical  compression  whereas 
zone  3  is  supposed  to  remain  unperturbed. 

less  marked  when  the  size  varies  from  10  to  15.  Interestingly  and 
contrary  to  the  capillary  pressure  curve,  the  size  effect  is  effective 
over  most  of  the  curve  and  not  only  in  the  vicinity  of  the  high  gas 
saturation.  Again  these  results  suggest  that  the  relative  permeabil¬ 
ities  should  be  determined  for  the  suitable  size  of  thin  system  and 
that  determination  over  systems  larger  than  the  one  actually  used 
in  the  application  (here  the  fuel  cell)  is  likely  to  introduce  significant 
errors. 

5.  Impact  of  reduced  exit  area  and  compression  effects 

As  sketched  in  Fig.  1,  only  a  fraction  0  of  the  GDL  unit  cell 
exit  surface,  i.e.  the  surface  of  the  GDL  in  contact  with  the  bi¬ 
polar  plate,  is  actually  in  contact  with  the  gas  channel.  Also,  due 
to  the  bi-polar  plate  grooved  structure,  it  is  expected  that  the  part 
of  the  GDL  in  contact  with  the  channels  is  less  mechanically  com¬ 
pressed  than  the  part  in  contact  with  the  solid  materials  when  the 
GDL  is  in  place  within  the  fuel  cell.  In  continuation  and  comple¬ 
ment  to  previous  studies,  e.g.  [9,21,22]  and  references  therein,  the 
possible  impact  of  these  effects  is  explored  in  this  section.  Before 
presenting  the  results  some  definitions  have  to  be  given.  Let  Lc  be 
the  fraction  of  L  corresponding  to  the  channel  (Fig.  1).  The  GDL 
unit  cell  exit  surface  open  fraction  is  then  given  by  0  =  Lc/L.  As 


O 


sketched  in  Fig.  7,  three  zones  are  distinguished  in  the  GDL  unit 
cell.  Zones  1  and  2  are  in  contact  with  the  bi-polar  plate  solid  mate¬ 
rials  whereas  zone  3  is  located  right  underneath  the  channel.  To 
simulate  the  effect  of  differential  compression,  the  method  is  as 
follows.  Starting  from  a  given  realization  of  the  GDL  unit  cell,  it 
is  assumed  that  the  compression  only  affects  zones  1  and  2.  In 
these  zones,  it  is  assumed  that  the  main  effect  of  compression  (at 
least  as  regards  the  impact  of  compression  on  the  liquid  invasion) 
is  to  reduce  the  size  d  of  a  throat  or  a  pore  according  to  the  sim¬ 
ple  rule:  dCOmp=n'd  where  a  is  the  compression  coefficient  with 
0<o'<  1.  Hence  a  =  \  corresponds  to  the  case  without  differential 
compression.  All  the  results  presented  in  this  section  are  obtained 
for  N=  10. 

5.1.  Overall  saturation  at  breakthrough  and  saturation  profiles 

Fig.  8a  shows  the  evolution  of  the  average  overall  saturation  in 
the  unit  cell  as  well  as  the  overall  saturation  in  the  three  zones 
as  a  function  of  0  at  breakthrough  for  a  =  1  (no  differential  com¬ 
pression).  As  expected  (5)Bt  increases  with  the  fraction  of  the  exit 
surface  in  contact  with  the  solid  materials.  The  saturation  tends 
to  be  slightly  higher  in  the  zones  located  underneath  the  solid 
materials  (zones  1  and  2).  The  effect  of  differential  compression 
is  depicted  in  Fig.  8b  for  0  =  0.5.  We  recall  that  the  invading  fluid 
invades  preferentially  the  interfacial  bonds  of  largest  size  (i.e.  of 
smallest  capillary  pressure  threshold)  in  the  invasion  percolation 
regime.  As  a  result,  the  compressed  zones  are  not  visited  at  all  with 
a  sufficient  compression.  This  suggests  that  the  differential  com¬ 
pression  can  be  beneficial:  liquid  water  would  be  transported  in 
zone  3  (the  less  compressed  zone)  whereas  the  most  compressed 
zones  (zones  1  and  2)  would  contain  no  or  significantly  less  water 
and  therefore  would  offer  zones  where  the  gas  transport  is  eas¬ 
ier. 

Fig.  9  shows  how  the  saturation  profile  (along  the  thickness)  is 
affected  by  the  reduced  exit  area  and  the  differential  compression 
effect.  These  profiles  are  obtained  by  computing  the  saturation  over 
successive  slices  of  network  for  many  realizations  (the  error  bars  in 
Fig.  9  correspond  to  ±1  standard  deviation  around  the  mean  values). 
The  profiles  in  Fig.  9  have  the  concave  shape  typical  of  IP  profiles  at 
breakthrough,  e.g.  [9]  and  references  therein. 

The  effect  of  differential  compression  shown  in  Figs.  8b  and  9b 
is  fully  consistent  with  what  is  expected  if  one  considers,  as  here, 
that  the  net  result  of  differential  compression  is  purely  geometri¬ 
cal  and  leads  to  more  narrow  throats  and  pores  in  the  compressed 
zones.  The  situation  may  be  actually  subtler  if  one  considers  that 
the  in-plane  compression  is  different  from  the  through  plane  com- 


Fig.  8.  (a)  Evolution  of  overall  saturation  at  breakthrouh  in  the  different  zones  as  a  function  of  free  outlet  surface  fraction  (without  differential  compression,  i.e.  a  =  1),  (b) 
evolution  of  overall  saturation  at  breakthrouh  as  a  function  of  differential  compression  factor  a  (for  a  free  outlet  surface  fraction  0  =  0.5). 
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Fig.  9.  Evolution  of  saturation  profile  at  breakthrough  (a)  Effect  of  free  outlet  surface  fraction  (a  =  1 ).  (b)  Effect  of  differential  compression  (0  =  0.5). 


pression  and  also  if  compression  effects  induced  modification  of 
local  wettability  conditions  as  suggested  in  [23]. 

5.2.  Capillary  pressure 

The  effect  of  differential  compression  on  the  apparent  capil¬ 
lary  pressure  is  shown  in  Fig.  10.  As  can  be  seen  from  Fig.  10,  the 
capillary  pressure  curve  is  not  affected  for  the  high  gas  saturation. 
The  capillary  pressure  increases  with  the  compression  in  the  range 
of  intermediate  saturations.  This  is  obviously  consistent  with  the 
fact  that  the  sizes  of  throats  located  in  the  compressed  zones  are 
assumed  to  become  narrower  with  the  compression.  Changing  the 
outlet  open  surface  fraction  0  in  the  absence  of  compression  has 
no  significant  effect  on  PC(S)  if  one  assumes,  as  we  do  in  this  paper, 
that  trapping  phenomena  can  be  neglected. 

5.3.  Permeability  and  relative  permeabilities 

Fig.  11  shows  the  impact  of  reduced  exit  area  and  compression 
effect  on  apparent  permeability.  Contrary  to  the  intrinsic  perme¬ 
ability  of  the  porous  medium,  which  is  supposed  to  be  a  local 
property  of  the  porous  microstructure,  the  apparent  permeabil¬ 
ity  characterizes  the  global  permeability  of  the  system  formed  by 
the  GDL  unit  cell  in  the  presence  of  the  bipolar  plate  and  associ¬ 


ated  effects  of  compression  and  reduced  exit  area.  The  apparent 
permeability  can  therefore  be  defined  as 


where  Q.  is  the  flow  rate,  Pen  the  pressure  at  the  entrance  of  the 
system  and  Pex  the  average  pressure  at  the  exit,  i.e.  in  the  channel 
at  z=t. 

As  discussed  in  [23],  the  apparent  permeability  decreases  when 
the  exit  area  is  reduced  owing  to  the  deformation  of  streamlines 
within  the  porous  domain.  This  is  exemplified  in  Fig.  11  a  for  the  par¬ 
ticular  system  considered  here.  As  depicted  in  Fig.  lib,  the  effect  of 
differential  compression  on  apparent  permeability  is  not  important 
over  the  range  of  compression  factor  investigated. 

As  shown  in  Fig.  12,  the  situation  is  quite  different  as  regards 
the  apparent  relative  permeabilities  (defined  as  fcra  =  {ftQjI<a{Pen  - 
Pex)){£/L2),  where  /x,  Q.  and  P  are  the  dynamic  viscosity,  flow  rate 
and  pressure  in  the  considered  fluid),  i.e.  the  main  factor  affecting 
the  apparent  relative  permeabilities  is  the  compression.  This  was 
expected  since  we  know  from  Section  5.1  that  the  differential  com¬ 
pression  has  a  very  significant  effect  on  the  liquid  water  distribution 
during  the  invasion  with  preferential  invasion  of  the  zone  located 
under  the  channel  for  a  sufficient  compression. 


Fig.  10.  Evolution  of  capillary  pressure  curve  with  compression  factor  a  (0  =  0.5). 


6.  Injection  from  an  isolated  source 

For  all  the  two-phase  flow  simulations  presented  so  far,  we 
have  assumed  that  liquid  water  was  present  over  the  entire  unit 
cell  entrance  surface  and  also  supposed  implicitly  that  the  liquid 
pressure  was  uniform  over  this  surface.  As  discussed  in  Section  3, 
this  case  can  be  representative  of  ex-situ  experimental  tests  aim¬ 
ing  at  characterizing  the  GDL  properties.  However,  the  boundary 
condition  for  a  GDL  in  a  fuel  cell  could  be  different  since  water 
has  to  travel  through  the  active  layer  and  the  MPL,  if  any,  before 
reaching  the  GDL.  Assume  for  example  the  presence  of  a  MPL.  The 
MPL  pores  are  much  smaller  than  GDL  ones.  Assume  for  simplicity 
that  the  water  transfer  within  the  MPL  can  also  be  described  by 
IP  (we  recall  that  the  MPL  is  hydrophobic  which  is  consistent  with 
the  IP  regime).  Then  water  should  break  at  the  MPL/GDL  interface 
through  one  pore  of  the  MPL.  According  to  this  extreme  case,  this 
would  imply  that  water  enters  the  GDL  unit  cell  through  only  one 
pore  of  the  GDL  entrance  surface.  It  is  likely  that  water  enters  in 
fact  through  several  isolated  pores  since  water  is  produced  from 
many  independent  sources  (^the  platinum  grains)  in  the  active 
layers,  but  the  details  are  still  an  open  question.  For  the  moment, 
we  simply  wish  to  explore  the  impact  of  a  single  injection  source. 
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Fig.  11.  (a)  Evolution  of  apparent  permeability  as  a  function  of  free  outlet  surface  fraction  (a  =  1 ).  (b)  Evolution  of  apparent  permeability  as  a  function  of  differential  compression 
(0  =  0.5). 


Hence  we  assume  that  water  invades  the  GDL  from  a  single  throat 
arbitrarily  chosen  as  the  throat  located  in  the  middle  of  the  unit 
cell  entrance  surface.  The  water  invasion  simulations  are  still  based 
on  the  IP  algorithm  and  have  been  performed  for  <2>  =  1  and  a  =  1, 
i.e.  no  restricted  exit  area,  no  compression.  As  shown  in  Fig.  4,  the 
injection  boundary  condition  has  a  significant  impact  on  the  overall 
saturation  at  breakthrough.  As  expected  a  smaller  fraction  of  the 
pore  space  is  invaded  when  water  enters  from  an  isolated  local 
source.  As  can  be  seen  from  Fig.  4,  (S)bt  increases  quasi-linearly 
with  the  unit  cell  thickness  N  for  the  one  point  injection  mode. 

The  saturation  profiles  at  breakthrough  for  the  two  injection 
modes  are  compared  in  Fig.  13.  Contrary  to  surface  injection,  the 
one-point  injection  leads  to  non-concave,  non-monotonic  profiles 
except  for  the  thinnest  systems  considered  for  which  the  profile  is 
almost  flat.  Hence  the  concave  shape  observed  in  the  traditional 
IP  simulations  is  strongly  dependent  on  the  boundary  condition 
and  this  can  contribute  to  explain  why  profiles  deduced  from 
in  situ  experiments,  e.g.  [24,25]  are  not  concave.  Also,  it  is  clear 
from  Fig.  13  that  the  one-point  injection  mode  is  significantly 
more  favourable  to  gas  transport  since  the  gas  saturation  is  higher 
everywhere  along  the  profile. 

7.  Continuum  model  vs  pore  network  model  steady  state 
saturation  profiles 

In  this  section  we  consider  again  the  basic  situation  (0  =  1  and 
a  =  1 )  and  compare  the  steady  state  saturation  profile  at  break¬ 


through  predicted  on  the  one  hand  thanks  to  the  IP  algorithm  on 
the  pore  network  and  on  the  other  hand  with  the  traditional  con¬ 
tinuum  model.  Following  an  approach  similar  to  the  one  presented 
in  [2],  the  continuum  model  leads  to 

IM  _  I<krl  dPc  dS 

2  F  v  dS  dz  1  J 

where  /<r£  is  the  water  relative  permeability.  Using  i  as  reference 
length,  Eq.  (8)  is  expressed  in  dimensionless  form  as 


dz  VF  Vs  kri(dJ/dS)  K  J 

where  J  is  the  Leverett  function  (Pc  =  Pwater  -  Pgas  - 
(ycosfl/v'K/eJKS)).  Eq.  (9)  is  solved  using  a  classical  fourth 
order  Runge-Kutta  method  in  conjunction  with  the  boundary 
condition  5  =  5*  in  z  =  l  (Z  =  1 ).  As  an  example,  we  consider  the  case 
of  the  40x40x10  system.  Numerical  fits  of  7(5)  and  kr^(S)  are 
deduced  from  the  data  reported  in  Figs.  5  and  6  for  this  system. 
One  problem  is  to  specify  S\  i.e.  the  boundary  condition  atZ=  1.  An 
obvious  option  is  to  take  S*  =  (5)Bt  since  kTi(S)  >  0  only  for  S  >  (S)B t- 
Numerical  tests  show,  however,  that  the  saturation  profile  that  is 
determined  from  Eq.  (9)  is  not  sensitive  to  the  exact  value  chosen 
for  5*  provided  that  the  stepsize  is  sufficiently  small  and  5*  is  close 
to  (S)bt. 

The  comparison  between  the  saturation  profiles  determined 
using  the  IP  algorithm  and  the  one  deduced  from  Eq.  (9)  (for  two 
values  of  capillary  number)  is  depicted  in  Fig.  14.  As  can  be  seen 


Fig.  12.  Evolution  of  apparent  relative  permeabilities  as  a  function  of  gas  overall  saturation  (a)  effect  of  free  outlet  surface  fraction  (a  =  1 ),  (b)  effect  of  differential  compression 
(0  =  0.5). 
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Fig.  13.  Influence  of  injection  boundary  condition  on  the  evolution  of  slice  average  liquid  saturation  along  the  network  thickness,  (a)  Surface  injection  and  (b)  one  point 
injection. 


from  Fig.  14,  there  are  very  significant  differences  between  the  IP 
profiles  and  the  ones  predicted  with  the  continuum  model.  Note 
in  particular  the  convex  shape  of  profiles  and  above  all  the  signif¬ 
icantly  higher  saturations  predicted  in  the  second  half  of  profiles 
with  the  continuum  model.  This  was  expected  since  it  is  widely 
admitted  that  the  continuum  model  cannot  describe  the  capillary 
fingering  regime  corresponding  to  IP.  It  is  however  interesting  to  see 
quantitatively  in  this  example  how  wrong  is  the  continuum  model. 
This  is  discussed  further  in  the  next  section. 

8.  Discussion 

As  illustrated  in  Fig.  14,  predicting  the  GDL  pore  blockage  by 
the  water  using  the  traditional  continuum  two-phase  flow  model 
is  likely  to  lead  to  serious  errors.  One  first  problem  is  that  the  water 
invasion  regime  in  a  GDL  is  expected  to  correspond  to  the  capil¬ 
lary  fingering  regime  (see  Section  2).  This  regime  leads  to  fractal 
invasion  patterns  and  cannot  be  described  using  the  traditional 
continuum  approach  (the  continuum  approach  to  porous  media  is 


0  0,2  0,4  0,6  0,8  1 


(Z-a)fafN 

Fig.  14.  Comparison  of  saturation  profiles  predicted  using  the  pore  network  model 
(PNM)  for  the  two  boundary  conditions  sketched  in  Fig.  1  and  the  profile  predicted 
with  the  continuum  model. 


conceptually  based  on  the  existence  of  Representative  Elementary 
Volume,  e.g.  [26],  and  such  a  REV  cannot  be  defined  with  a  fractal 
distribution  of  the  fluids). 

A  related  problem  is  that  a  GDL  is  very  thin  (measured  in  pore 
size).  As  a  result,  it  becomes  questionable  to  use  the  traditional 
continuum  model  as  a  local  model.  For  example  the  parameters 
s,  K,  kr£,  ( dJ/dS )  in  Eq.  (8)  are  traditionally  considered  as  local 
parameters,  i.e.  defined  over  a  REV.  In  the  case  of  the  GDL  unit 
cell  considered  in  this  paper,  we  face  the  odd  situation  in  which 
the  system  size  is  on  the  order  or  even  lower  than  the  REV  size. 
Under  these  circumstances,  it  becomes  clearly  doubtful  to  give 
its  usual  local  meaning  to  the  classical  two-phase  flow  model. 
Hence,  to  obtain  the  profiles  shown  in  Fig.  14,  we  have  solved 
Eq.  (8)  as  a  local  equation  using  parameters  determined  over  the 
entire  system,  which  is  not  at  all  in  agreement  with  the  tradi¬ 
tional  concept  of  length  scale  separation  underlying  the  continuum 
model.  The  scale  dependence  of  the  parameters  exemplified  in 
Section  4  is  another  illustration  of  the  lack  of  REV  in  the  present 
context. 

The  fundamental  differences  between  the  continuum  models 
and  the  IP  approach  can  also  be  seen  from  the  fact  that  the  contin¬ 
uum  model  profiles  are  capillary  number  dependent  whereas  there 
are  independent  of  viscous  effects  with  the  IP  model  (see  Section 
2).  Also,  as  pointed  out  in  [17],  a  boundary  condition  should  be 
specified  at  the  GDL/channel,  i.e.  the  “local”  saturation,  with  the 
continuum  model  whereas  the  IP  model  predicts  this  saturation, 
which  therefore  need  not  to  be  specified  (this  implies,  however, 
that  there  is  no  limitation  in  the  channel  for  evacuating  the  water 
reaching  the  channel  from  the  GDL). 

Owing  to  the  lack  of  length  scale  separation,  using  the  con¬ 
tinuum  model  as  a  local  model  in  2D  or  3D  simulations  so  as  to 
investigate  for  example  the  effect  of  the  outlet  restricted  area  at 
the  GDL/channel  interface  is  likely  to  lead  to  predictions  still  poorer 
than  for  the  ID  case  considered  in  Section  7. 

Owing  to  the  poor  performances  of  continuum  models,  pore  net¬ 
work  models  offer  an  attractive  alternative.  Yet  there  are  several 
problems.  As  illustrated  in  Figs.  13  and  14,  the  IP-PNM  prediction  is 
quite  sensitive  to  the  injection  boundary  condition.  Also,  the  con¬ 
ditions  used  in  this  paper  lead  to  a  single  breakthrough  point  at  the 
GDL/channel  interface,  which  is  not  in  agreement  with  the  visuali¬ 
sation  studies,  e.g.  [27],  which  rather  indicate  several  breakthrough 
points,  i.e.  the  formation  of  several  droplets  and  not  a  single  one 
over  channel  distances  comparable  to  the  GDL  unit  cell  size  con¬ 
sidered  in  the  present  study.  If  viscous  effects  can  be  discarded  as 
discussed  in  Section  2,  this  can  be  explained  by  either  condensa¬ 
tion  effects  within  the  GDL  since  the  GDL  is  cooler  than  the  active 
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layer  or  by  the  fact  that  liquid  water  enters  the  GDL  through  several 
independent  injection  points,  [9,17]. 

Owing  to  the  IP-PNM  model  sensitivity  to  the  injection  boundary 
condition  (see  Figs.  13  and  14),  the  question  of  the  proper  injection 
condition  to  be  imposed  at  the  GDL  inlet  surface  appears  as  very 
important.  As  outlined  in  [17],  it  seems  difficult  to  specify  this  con¬ 
dition  without  the  full  consideration  of  the  water  production  and 
transfers  in  the  GDL  adjacent  layers,  i.e.  the  active  layer  and  the 
MPL.  Also,  it  is  widely  admitted  that  an  increase  in  the  current  den¬ 
sity  should  lead  to  an  increase  in  the  overall  saturation  of  the  GDL. 
Within  the  IP-PNM  framework,  this  would  mean  an  increase  in  the 
number  of  injection  points. 

Throughout  this  study  we  have  assumed  that  the  GDL  was 
hydrophobic  and  that  invasion  could  be  simulated  using  the  IP  algo¬ 
rithm.  As  pointed  out  in  [4],  the  IP  algorithm  is  strictly  valid  only 
for  sufficiently  high  contact  angle  (measured  in  the  invading  phase). 
For  instance,  the  simulations  shown  in  [4]  for  a  2D  model  system 
indicates  that  the  invasion  pattern  is  not  strictly  an  IP  one  for  con¬ 
tact  angles  representative  of  water  on  PTFE  owing  to  changes  in 
the  liquid-gas  interface  local  growth  mechanisms,  see  [4]  for  more 
details.  The  range  of  contact  angles  for  which  the  IP  model  becomes 
in  serious  error  can  be  perfectly  characterized  for  systems  similar 
as  the  2D  ones  considered  in  [4].  For  a  much  more  complex  3D 
microstructure  like  a  GDL,  this  is  an  unexplored  question. 

Also,  it  seems  to  be  more  correct  to  consider  a  GDL  as  a  system 
of  mixed  wettability  rather  than  a  purely  hydrophobic  system,  e.g. 
[7,28,29].  Liquid  invasion  dominated  by  capillary  effects  in  a  system 
of  mixed  wettability  cannot  be  described  using  the  IP  algorithm 
and  more  complex  local  invasion  rules  must  be  used,  e.g.  [4,7,29]. 
Flence,  this  aspect  should  be  taken  into  account  in  the  pore  net¬ 
work  simulations.  However,  it  is,  by  no  means,  expected  that  the 
consideration  of  mixed  wettability  effects  will  change  the  conclu¬ 
sions  of  the  present  study  regarding  the  poor  predictions  that  can 
be  expected  from  continuum  models.  The  invasion  process  will  be 
still  dominated  by  capillary  effects  and  the  lack  of  length  scale  sep¬ 
aration  will  be  still  present.  We  have  also  ignored  the  non-isotropic 
nature  of  the  GDL  microstructure.  However,  it  is  clear  that  almost  all 
GDL  types  are  highly  anisotropic.  Although  this  will  not  change  the 
main  results  of  the  present  study,  the  influence  of  GDL  anisotropy 
on  liquid  invasion  would  deserve  to  be  studied  in  some  details  in  a 
future  work. 

9.  Conclusion 

As  illustrated  in  this  paper,  the  traditional  continuum  mod¬ 
els  widely  used  for  describing  two-phase  flow  in  GDLs  should  be 
considered  as  leading  to  only  poor  approximations  of  water  distri¬ 
bution  in  the  GDL. 

The  main  reasons  of  the  continuum  model  poor  performances 
are  the  lack  of  length  scale  separation  (which,  among  other  things, 
leads  to  scale  dependent  transport  properties)  and  the  fact  that  the 
invasion  regime  is  dominated  by  capillary  effects  (this  introduces 
a  second  lack  of  length  scale  separation  owing  to  the  phases  frac¬ 
tal  distribution  corresponding  to  the  capillary  fingering  regime).  In 
other  terms,  the  condition  underlying  the  continuum  approach  to 
porous  media,  i.e.  the  existence  of  a  representative  elementary  vol¬ 


ume  of  size  much  smaller  than  the  porous  domain,  is  not  met  in  a 
GDL. 

Due  to  the  limitations  of  continuum  models,  pore  network 
models  (PNMs)  appear  as  a  good  trade-off  between  accuracy  and 
computational  effort. 

However,  several  aspects  related  to  the  use  of  PNMs  needs  to  be 
clarified.  This  notably  includes  the  injection  condition  and  a  better 
assessment  of  wettability  related  effects  (relevance  of  IP  model  for 
contact  angles  representative  of  water  on  PTFE,  mixed  wettability). 

Also,  since  the  invasion  pattern  is  strongly  dependent  on  the 
injection  boundary  condition,  this  study  suggests  that  a  proper 
analysis  of  water  transport  in  a  GDL  cannot  be  performed  with¬ 
out  consideration  of  water  generation  and  transport  in  the  adjacent 
porous  layers  (active  layer  and  MPL). 
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